! 这是C3D8(I)单元的主程序文件
program sdyna
    use iso_fortran_env, only: real64
    use mesh_reader_mod
    use dbc_reader_mod
    use fbc_reader_mod
    use amplitude_reader_mod
    use writer_mod
    use solid_mod,only: c3d8_ke,c3d8_me,c3d8i_ke,c3d8i_me,nndofs,get_ndof_inds
    use ls_mod,only: do_hex_static_solve,do_hex_eigs_solve8,do_hex_newmark
    use sconstants_mod
    implicit none
    ! 求解参数变量
    logical                             :: opened
    integer                             :: nnum, dbc_num, fbc_num , amps_num
    integer                             :: nev
    integer                             :: i,j,iostat, fileunit
    integer                             :: gdof, eltype_code, ninds(nndofs)
    character(len=2)                    :: solve_type
    real(real64)                        :: alpha,beta  ! 瑞利阻尼系数
    real(real64)                        :: inc_time,step_time,E,nu,rho
    ! 输入数据/文件变量
    character(len=256)                  :: nodefp,elemfp,dbcfp,fbcfp,ampfp,configfp,inpfile,elemtype
    integer,allocatable                 :: el_conn(:,:),dbc_pos(:,:),dbc_amp(:),fbc_pos(:,:),fbc_amp(:)
    real(real64),allocatable            :: nodes(:,:),dbc_val(:),fbc_val(:), amps_array(:,:,:)
    ! 静力分析结果变量
    real(real64),allocatable            :: nodal_u(:,:)
    real(real64),allocatable            :: nodal_f(:,:)
    ! 模态分析结果变量
    real(real64),allocatable            :: eigvals(:)
    real(real64)                        :: eigfreq
    real(real64),allocatable            :: eigvecs(:,:)

    amps_num    = 0
    inc_time    = 0.0_real64
    step_time   = 0.0_real64
    E           = 0.0_real64
    nu          = 0.0_real64
    rho         = 0.0_real64
    alpha       = 0.0_real64
    beta        = 0.0_real64
    nev         = 0
    eigfreq     = 0.0_real64
    ! ==========================================================
    !                       read inputs file
    ! ==========================================================
    inquire(unit=100,opened=opened)
    if (.not. opened) then
        call get_command_argument(number=1,value=inpfile) ! 获取命令行参数
        open(unit=100,status='old',file=inpfile,action='read',iostat=iostat)
        if (iostat /=0) then
            print *, 'open commandline args file failed'
            stop
        end if
        read(100,'(A)') nodefp      ! 读取节点坐标文件名
        read(100,'(A)') elemfp      ! 读取单元连接文件名
        read(100,'(A)') dbcfp       ! 读取dbc文件名
        read(100,'(A)') fbcfp       ! 读取fbc文件名
        read(100,'(A)') ampfp       ! 读取amp文件名
        read(100,'(A)') configfp    ! 读取配置文件名 
        close(100) ! 关闭文件
    end if

    ! 读取 config.txt 数据
    print *, '==> Start reading confg file'
    open(newunit=fileunit,file=configfp,status='old',action='read',iostat=iostat)
    read(fileunit,*) solve_type,elemtype                        ! 读取分析步类型,单元类型

    print *, '+ --> solve type:',solve_type
    print *, '+ --> element type:',elemtype
    elemtype = trim(adjustl(elemtype))                          ! 去除空格
    if (elemtype=='c3d8') then
        eltype_code = 1
    else if ( elemtype=='c3d8i' ) then
        eltype_code = 2
    end if
    
    select case(solve_type)
    case(linear_static_step)
        print *, '==> Start reading input data !'
        call load_nodes(nodefp,nodes,nnum,iostat)           ! 读取 node.txt 数据
        if (iostat > 0) print *, "load node failed, iostat=", iostat
        call load_elements(elemfp,8,el_conn,iostat)         ! 读取 element.txt 数据
        if (iostat > 0) print *, "load element failed, iostat=", iostat
        call load_dbcs_ls(dbcfp,dbc_pos,dbc_val,dbc_num,iostat)! 读取 dbc.txt 数据
        if (iostat > 0) print *, "load dbc failed, iostat=", iostat
        call load_fbcs_ls(fbcfp,fbc_pos,fbc_val,fbc_num,iostat) ! 读取 fbc.txt 数据
        if (iostat > 0) print *, "load fbc failed, iostat=", iostat

        read(fileunit,*) E, nu                              ! 读取材料参数
        close(fileunit)                                     ! 关闭配置文件
        print *, "==>Start linear static analysis"
        allocate(nodal_u(nnum,nndofs))                      ! 结果数组分配内存
        allocate(nodal_f(nnum,nndofs))                      ! 结果数组分配内存
        ! 调用线性分析子程序
        gdof = nnum*nndofs                                  ! 全局自由度
        call do_hex_static_solve(nodes, el_conn, e, nu, eltype_code,&
                                dbc_pos, dbc_val, fbc_pos, fbc_val,&
                                gdof, nodal_u, nodal_f)
        ! 写入tecplot结果文件
        call get_command_argument(number=2, value=inpfile)  ! 获取结果文件名
        open(newunit=fileunit, file=inpfile, status='replace',action='write',iostat=iostat)
        if (iostat/=0) print *, 'open commandline args file failed'
        write(fileunit,'(A)') 'TITLE = "Linear Static Analysis Result"'
        write(fileunit,'(A)') 'VARIABLES = "X","Y","Z","U1","U2","U3","RF1","RF2","RF3"'
        write(fileunit,'(A,I8,A,I10,A)') 'ZONE T="SOLID_MESH",N=',nnum,',E=',size(el_conn,1),',ZONETYPE=FEBrick,DATAPACKING=POINT'
        do i=1,nnum                                         ! 写入节点和节点数据
            write(fileunit,'(3F28.10,3F30.16,3F30.16)') nodes(i,1),nodes(i,2),nodes(i,3),&
            nodal_u(i,1),nodal_u(i,2),nodal_u(i,3),nodal_f(i,1),nodal_f(i,2),nodal_f(i,3)
        end do
        
        do i=1,size(el_conn,1)                              ! 写入单元连接定义
            write(fileunit,'(8I16)') el_conn(i,:)
        end do
        
        close(fileunit)
        print *, "==>End linear static step"
        print *, "==>Results are saved to ", inpfile
        deallocate(nodal_u,nodal_f)                         ! 释放结果数组内存
        deallocate(nodes,dbc_val,fbc_val)                   ! 释放节点/位移和力约束的内存
        deallocate(el_conn,dbc_pos,fbc_pos)
    case(natural_frequency_step)
        print *, '==> Start reading input data !'
        call load_nodes(nodefp,nodes,nnum,iostat)           ! 读取 node.txt 数据
        if (iostat > 0) print *, "load node failed, iostat=", iostat
        call load_elements(elemfp,8,el_conn,iostat)         ! 读取 element.txt 数据
        if (iostat > 0) print *, "load element failed, iostat=", iostat
        call load_dbcs_nf(dbcfp,dbc_pos,dbc_num,iostat)! 读取 dbc.txt 数据
        if (iostat > 0) print *, "load dbc failed, iostat=", iostat

        read(fileunit,*) nev                                ! 读取特征值数
        read(fileunit,*) E, nu, rho                         ! 读取材料参数
        close(fileunit)                                     ! 关闭配置文件
        print *, "==>Start Natrual Frequency analysis"
        gdof = nnum * nndofs                                ! 全局节点自由度
        allocate(eigvals(nev))                              ! 分配数组内存
        allocate(eigvecs(gdof,nev))                         ! 分配数组内存
        ! 调用特征值求解子程序
        call do_hex_eigs_solve8(nodes,el_conn,e,nu,rho,nev,eltype_code,&
                    dbc_pos,gdof,eigvals,eigvecs)
        do i=1,nev
            print *, 'mode -',i,' eigen value =',eigvals(i)
        end do 
        ! 写入到tecplot文件
        call get_command_argument(number=2, value=inpfile)  ! 获取结果文件名
        open(newunit=fileunit, file=inpfile, status='replace',action='write',iostat=iostat)
        if (iostat/=0) print *, 'open commandline args file failed'
        write(fileunit,*) 'TITLE = "Natural Frequency modalshape"'
        write(fileunit,'(A)') 'VARIABLES = "X","Y","Z","U1","U2","U3"'
        do i=1,nev                                          ! 循环写入模态数据
            eigfreq=sqrt(eigvals(i))/(2.0_real64*pi)
            write(fileunit,'(A,I4,A,I8,A,I10,A,F15.4)') 'ZONE T="Mode-',i,'",N=',nnum,',E=',size(el_conn,1),&
            ',ZONETYPE=FEBrick,DATAPACKING=POINT, SOLUTIONTIME=',eigfreq
            do j=1,nnum
                ninds=get_ndof_inds(j)
                write(fileunit,'(F28.10,F28.10,F28.10,ES25.13,ES25.13,ES25.13)') nodes(j,1),&
                nodes(j,2),nodes(j,3),eigvecs(ninds(1),i),eigvecs(ninds(2),i),eigvecs(ninds(3),i)
            end do
            do j=1,size(el_conn,1)                          ! 写入单元连接定义
                write(fileunit,'(8I16)') el_conn(j,:)
            end do
        end do
        print *, "==>End Natrual Frequency step"
        print *, "==>Results are saved to ", inpfile
        close(fileunit)
        deallocate(eigvals,eigvecs,nodes,el_conn,dbc_pos,dbc_val)
    case(transient_dynamics_step)
        print *, '==> Start reading input data !'
        call load_nodes(nodefp,nodes,nnum,iostat)                           ! 读取 node.txt 数据
        if (iostat > 0) print *, "load node failed, iostat=", iostat
        call load_elements(elemfp,8,el_conn,iostat)                         ! 读取 element.txt 数据
        if (iostat > 0) print *, "load element failed, iostat=", iostat
        call load_dbcs_dy(dbcfp,dbc_pos,dbc_val,dbc_amp,dbc_num,iostat)     ! 读取 dbc.txt 数据
        if (iostat > 0) print *, "load dbc failed, iostat=", iostat
        call load_fbcs_dy(fbcfp,fbc_pos,fbc_val,fbc_amp,fbc_num,iostat)     ! 读取 fbc.txt 数据
        if (iostat > 0) print *, "load fbc failed, iostat=", iostat
        call load_amplitude_dy(ampfp,amps_array,amps_num,iostat)                     ! 读取 amplitude.txt 数据
        if (iostat > 0) print *, "load amplitude failed, iostat=", iostat

        read(fileunit,*) step_time,inc_time                 ! 总时间, 时间增量
        read(fileunit,*) E, nu, rho                         ! 读取材料参数
        read(fileunit,*) alpha, beta                        ! 瑞利阻尼的alpha, beta
        close(fileunit)
        print*, '==> Reading material paramaters OK !'
        gdof=nnum*nndofs                                    ! 总自由度
        call get_command_argument(number=2, value=inpfile)  ! 获取结果文件名
        call do_hex_newmark(nodes,el_conn,e,nu,rho,alpha,beta,eltype_code,dbc_pos,&
        dbc_val,dbc_amp,fbc_pos,fbc_val,fbc_amp,amps_array,gdof,step_time,inc_time,inpfile)
        print *, "==>End Dynamic Analysis step"
        print *, "==>Results are saved to ", inpfile
        
        deallocate(amps_array,nodes,el_conn,dbc_pos,dbc_val,dbc_amp,fbc_pos,fbc_val,fbc_amp)
    end select
end program sdyna
